# ===================================================
# R codes for Figure 4-Figure 8
# ===================================================
library(dplyr)
library(tidyverse)
library(igraph)
library(ggraph)
library(qgraph)
library(ggthemes)
library(ggnetwork)
library(sna)
library(network)
library(GGally)
library(intergraph)
library(ggplot2)
library(dplyr)
library(lubridate)
library(car)
library(stargazer)
library(gridExtra)
library(tidyverse)
library(survival)
library(effects)
library(ggeffects)
library(interactions)
library(margins)
library(visreg)
library(tseries)


## ========================================================
## Figure 4: newly signed SPs
## ========================================================

df <- read_csv("ally sp count.csv")
library(readr)

# Select relevant columns and reshape to long format
annual_data <- df %>%
  select(year, country_name, ally_times, non_ally_times) %>%
  pivot_longer(cols = c(ally_times, non_ally_times),
               names_to = "type",
               values_to = "count") %>%
  mutate(type = dplyr::recode(type,
                              ally_times = "Ally",
                              non_ally_times = "Non-Ally"),
         type = factor(type, levels = c("Ally", "Non-Ally")))

# Create the faceted bar plot with custom black and light grey colors
figure4 <- ggplot(annual_data, aes(x = year, y = count, fill = type)) +
  geom_bar(stat = "identity", position = "dodge", color = "black", linewidth = 0.2) +
  facet_wrap(~ country_name, scales = "free_y") +
  scale_fill_manual(values = c("Ally" = "black", "Non-Ally" = "lightgrey")) +
  labs(x = "Year",
       y = "Strategic Partnerships Newly Signed (counts)",
       fill = NULL) +
  theme_classic(base_family = "Times", base_size = 12) +
  theme(
    strip.text = element_text(face = "bold", size = 12),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 11),
    legend.position = "bottom",
    legend.key.size = unit(0.4, "cm"),
    panel.spacing = unit(1, "lines"),
    plot.title = element_blank()
  )

figure4


## =============================================================
## Data preparation for Figure 5-Figure 8
## =============================================================

# Import the Data Set
att <- read_csv("global_sp_data.csv")

# Create lagged variables
att <- arrange(att, country, year)
att <- att %>%
  group_by(country) %>%
  mutate(prev_year_gdp = lag(as.numeric(GDP))) %>%
  mutate(prev_year_population = lag(as.numeric(population)))%>%
  mutate(prev_year_miexp = lag(as.numeric(`military expenditure`)))%>%
  mutate(prev_year_import = lag(as.numeric(import)))%>%
  mutate(prev_year_export = lag(as.numeric(export)))%>%
  mutate(prev_year_tradevol = lag(as.numeric(trade_volumn)))%>%
  mutate(prev_year_export.to.usa = lag(as.numeric(export_to_usa)))%>%
  mutate(prev_year_import.from.usa = lag(as.numeric(import_from_usa)))%>%
  mutate(prev_year_export.to.china = lag(as.numeric(export_to_china)))%>%
  mutate(prev_year_import.from.china = lag(as.numeric(import_from_china)))%>%
  mutate(prev_year_arms = lag(as.numeric(US_export_arms)))%>%
  mutate(prev_year_pr = lag(as.numeric(pr)))%>%    #political rights
  mutate(prev_year_cl = lag(as.numeric(cl)))%>%    # civil liberty
  mutate(prev_year_status = lag(as.numeric(status)))%>% #freedom status (freedom house)
  mutate(prev_year_fh.total.reversed = lag(as.numeric(fh_total_reversed)))%>% #the higher 
  ungroup()

# Create a binary dependent variable = signing SP with U.S. that year
att <- att %>% mutate(bi.sp.us= ifelse(usa_sp_signed==TRUE, 1,0))

att.2 <- att %>%
  group_by(country) %>%  # Group by country to ensure the lead function works within each country
  mutate(gdp_next_year = lead(GDP, 1)) 

# Convert 'US_ally' to a factor with specified labels
att.2$US_ally <- factor(att.2$US_ally, levels = c(0, 1), labels = c("Non-Ally", "Ally"))

att.2 <- att %>%
  group_by(country) %>%  # Group by country to ensure the lead function works within each country
  mutate(gdp_next_year = lead(GDP, 1),
         military_exp_next_year = lead(`military expenditure`, 1),
         trade_next_year = lead(trade_volumn, 1),
         export_to_usa_next_year = lead(export_to_usa, 1),
         import_from_usa_next_year = lead(import_from_usa,1),
         export_to_china_next_year = lead(export_to_china,1),
         import_from_china_next_year = lead(import_from_china,1),
         arms_next_year = lead(US_export_arms,1),
         arms_next_year_fixed = lead(US_export_arms_fixed,1),
         fh_total_reversed_next_year = lead(fh_total_reversed, 1))  # Shift the GDP to the next year

att.2$bi.sp.us <- as.factor(att.2$bi.sp.us)
att.2$US_ally <- as.factor(att.2$US_ally)

# Convert 'US_ally' to a factor with specified labels
att.2$US_ally <- factor(att.2$US_ally, levels = c(0, 1), labels = c("Non-Ally", "Ally"))

# ========================================================================================================
# IV = Cumulative SPs; DV= Arms Transfer
# ========================================================================================================

att.2$country <- as.factor(as.character(att.2$country))
att.2$year <- as.factor(as.character(att.2$year))

## TWFE 
model.6 <- lm(log(arms_next_year_fixed+1)~usa_sp_times*US_ally+population + GDP+country+year, data=att.2)
summary(model.6)
stargazer(model.6, type = "text")

## F Country
model.6.f.country <- lm(log(arms_next_year_fixed+1)~usa_sp_times*US_ally+population+ GDP+country, data=att.2)
summary(model.6.f.country)

## F Year
model.6.f.year <- lm(log(arms_next_year_fixed+1)~usa_sp_times*US_ally+population+ GDP+year, data=att.2)
summary(model.6.f.year)

## F Nothing
model.6.f.no <- lm(log(arms_next_year_fixed+1)~usa_sp_times*US_ally+(population)+ (GDP), data=att.2)
summary(model.6.f.no)

## Table A1: Regression Table for Arms Transfer
stargazer(model.6, model.6.f.country, model.6.f.year, model.6.f.no,  type = "text")


# =====================================================
# Figure 5: Marginal Effect on Arms Transfer with the U.S.
# =====================================================

library(effects)
library(lattice)

## TWFE 
model.6 <- lm(log(arms_next_year_fixed+1)~usa_sp_times*US_ally+population + GDP+country+year, data=att.2)
summary(model.6)
stargazer(model.6, type = "text")

# Step 1: Create the interaction effect
interaction_effect <- Effect(c("usa_sp_times", "US_ally"), model.6)

# Step 2: Rename factor levels
levels(interaction_effect$x$US_ally) <- c("Non-Ally", "Ally")

# Step 3: Extract lattice object
p <- plot(interaction_effect, 
          main = " ",
          xlab = "# of Strategic Partnership with the U.S.",
          ylab = "Predicted Log of U.S. Arms Transfer Next Year",
          given.values = NULL,
          multiline = FALSE,
          rug = TRUE,
          ci.style = "bands",
          plot = FALSE  
)

# Step 4: Define custom strip function
custom_strip <- function(which.given, which.panel, var.name, factor.levels, ...) {
  strip.default(which.given, which.panel,
                var.name = "",  # removes "US_ally ="
                factor.levels = factor.levels, ...)
}

# Step 5: Update the lattice object with custom strip
figure5 <- update(p, 
            strip = custom_strip,
            scales = list(
              x = list(relation = "free", alternating = FALSE),  # only bottom ticks
              y = list(relation = "same", limits = c(0.5, 3.6))  # same y scale across panels
            ),
            par.settings = list(
              axis.line = list(col = "black")
            ),
            between = list(x = 1))  # spacing between panels
print(figure5)


# =====================================================
# Figure 6: Marginal Effect on Trade with the U.S.
# =====================================================

## TWFE 
model.8 <- lm(log(trade_next_year+1)~usa_sp_times*US_ally+population + GDP+country+year, data=att.2)
stargazer(model.8, type = "text")
# Step 1: Create the interaction effect
interaction_effect <- Effect(c("usa_sp_times", "US_ally"), model.8)

# Step 2: Rename factor levels
levels(interaction_effect$x$US_ally) <- c("Non-Ally", "Ally")

# Step 3: Extract lattice object
p <- plot(interaction_effect, 
          main = " ",
          xlab = "# of Strategic Partnership with the U.S.",
          ylab = "Predicted Log of U.S. Trade Next Year",
          given.values = NULL,
          multiline = FALSE,
          rug = TRUE,
          ci.style = "bands",
          plot = FALSE  # <---- Don't plot yet, just return the object
)

# Step 4: Define custom strip function
custom_strip <- function(which.given, which.panel, var.name, factor.levels, ...) {
  strip.default(which.given, which.panel,
                var.name = "",  # removes "US_ally ="
                factor.levels = factor.levels, ...)
}

# Step 5: Update the lattice object with custom strip
figure6 <- update(p, 
            strip = custom_strip,
            scales = list(
              x = list(relation = "free", alternating = FALSE),  # only bottom ticks
              y = list(relation = "same", limits = c(23.4, 24.7))  # same y scale across panels
            ),
            par.settings = list(
              axis.line = list(col = "black")
            ),
            between = list(x = 1))  # spacing between panels
print(figure6)


# =====================================================
# Figure 7: Marginal Effect on Import from the U.S.
# =====================================================

## TWFE 
model.9 <- lm(log(import_from_usa_next_year+1)~usa_sp_times*US_ally+population + GDP+country+year, data=att.2)
stargazer(model.9, type = "text")

# Step 1: Create the interaction effect
interaction_effect <- Effect(c("usa_sp_times", "US_ally"), model.9)

# Step 2: Rename factor levels
levels(interaction_effect$x$US_ally) <- c("Non-Ally", "Ally")

# Step 3: Extract lattice object
p <- plot(interaction_effect, 
          main = " ",
          xlab = "# of Strategic Partnership with the U.S.",
          ylab = "Predicted Log of Import from the U.S. Next Year",
          given.values = NULL,
          multiline = FALSE,
          rug = TRUE,
          ci.style = "bands",
          plot = FALSE  # <---- Don't plot yet, just return the object
)

# Step 4: Define custom strip function
custom_strip <- function(which.given, which.panel, var.name, factor.levels, ...) {
  strip.default(which.given, which.panel,
                var.name = "",  # removes "US_ally ="
                factor.levels = factor.levels, ...)
}

# Step 5: Update the lattice object with custom strip
figure7 <- update(p, 
            strip = custom_strip,
            scales = list(
              x = list(relation = "free", alternating = FALSE),  # only bottom ticks
              y = list(relation = "same", limits = c(11.6, 13.6))  # same y scale across panels
            ),
            par.settings = list(
              axis.line = list(col = "black")
            ),
            between = list(x = 1))  # spacing between panels
print(figure7)

# =====================================================
# Figure 8: Marginal Effect on Export to the U.S.
# =====================================================
## TWFE 
model.10 <- lm(log(export_to_usa_next_year+1)~usa_sp_times*US_ally+population + GDP+country+year, data=att.2)
summary(model.10)

# Step 1: Create the interaction effect
interaction_effect <- Effect(c("usa_sp_times", "US_ally"), model.10)

# Step 2: Rename factor levels
levels(interaction_effect$x$US_ally) <- c("Non-Ally", "Ally")

# Step 3: Extract lattice object

p <- plot(interaction_effect, 
          main = " ",
          xlab = "# of Strategic Partnership with the U.S.",
          ylab = "Predicted Log of Export to the U.S. Next Year",
          given.values = NULL,
          multiline = FALSE,
          rug = TRUE,
          ci.style = "bands",
          plot = FALSE)

# Step 4: Define custom strip function
custom_strip <- function(which.given, which.panel, var.name, factor.levels, ...) {
  strip.default(which.given, which.panel,
                var.name = "",  # removes "US_ally ="
                factor.levels = factor.levels, ...)
}

p <- update(p, 
            strip = custom_strip,
            scales = list(
              x = list(relation = "free", alternating = FALSE),  # only bottom ticks
              y = list(relation = "same", limits = c(11.8, 13.2))  # same y scale across panels
            ),
            par.settings = list(
              axis.line = list(col = "black")
            ),
            between = list(x = 1))  # spacing between panels

# Step 5: Update the lattice object with custom strip
figure8 <- update(p, strip = custom_strip)

print(figure8)


## ================================================
## DCA  
## ================================================
dca <- read_csv("04_dca.csv")
dca <- dca %>% rename(iso2c =`iso2c (A)`)
dca <- dca %>% mutate(country=`Country Name (A)`)
dca_usa <- dca %>% filter(`Country B`=="USA")                   
             
# Select the columns you want to merge from dca_usa
dca_usa_subset <- dca_usa[, c("iso2c", "country", "year", 
                              "dcaGeneralV1", "dcaGeneralV2", 
                              "dcaSectorV1", "dcaSectorV2", 
                              "dcaAnyV1", "dcaAnyV2", 
                              "terminated_ynever", "terminated_y10", 
                              "terminated_y5", "terminated_y3")]

# Now merge into att.2
att.2_merged <- merge(att.2, dca_usa_subset, 
                      by = c("country", "year"), 
                      all.x = TRUE)  # Keep all rows in att.2                   
att.2_merged$year <- as.numeric(as.character(att.2_merged$year))                   

att.2_merged <- att.2_merged %>% filter(year<=2010)                   

# Model: DV=arms/IV=SP & DCA

## dcaAnyV1
model.arm.sp.dcav1 <- lm(log(arms_next_year_fixed+1)~usa_sp_times+dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
summary(model.arm.sp.dcav1)
stargazer(model.arm.sp.dcav1, type = "text") # sp=0.461***/ dcaAnyV1= 1.699***

## sp*dcaAnyV1
model.arm.sp.dcav1.interaction <- lm(log(arms_next_year_fixed+1)~usa_sp_times*dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
summary(model.arm.sp.dcav1.interaction)
stargazer(model.arm.sp.dcav1.interaction, type = "text") # sp = 0.792***/ dcaAnyV1=1.892***/ sp*dcaAnyV1 = 0.801**

## dcaAnyV2
model.arm.sp.dcav2 <- lm(log(arms_next_year_fixed+1)~usa_sp_times+dcaAnyV2+US_ally+population + GDP+factor(year), data=att.2_merged)
summary(model.arm.sp.dcav2) # sp=0.463***/ dcaAnyV2=1.165***
stargazer(model.arm.sp.dcav1, model.arm.sp.dcav2, type = "text") 

## sp*dcaAnyV2
model.arm.sp.dcav2.interaction <- lm(log(arms_next_year_fixed+1)~usa_sp_times*dcaAnyV2+US_ally+population + GDP+factor(year), data=att.2_merged)
summary(model.arm.sp.dcav2.interaction) # sp=1.091***/ dcaAnyV2=1.349***/ sp*dcaAnyV2= -0.925***
stargazer(model.arm.sp.dcav1.interaction, model.arm.sp.dcav2.interaction, type = "text")   


# Model: DV=trade/IV=SP & DCA                   

## dcaAnyV1
model.trade.sp.dcav1 <- lm(log(trade_next_year+1)~usa_sp_times+dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
summary(model.trade.sp.dcav1)
stargazer(model.trade.sp.dcav1, type = "text") # sp=0.345***/ dcaAnyV1= 1.619***

## sp*dcaAnyV1
model.trade.sp.dcav1.interaction <- lm(log(trade_next_year+1)~usa_sp_times*dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
summary(model.trade.sp.dcav1.interaction)
stargazer(model.trade.sp.dcav1.interaction, type = "text") # sp = 0.468***/ dcaAnyV1=1.692***/ sp*dcaAnyV1 = -0.242**


# Model: DV=arms transfer / IV=terminated_y3,5,10

model.arm.sp.y3 <- lm(log(arms_next_year_fixed+1)~terminated_y3+dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
stargazer(model.arm.sp.y3, type="text")  

model.arm.sp.y5 <- lm(log(arms_next_year_fixed+1)~terminated_y5+dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
stargazer(model.arm.sp.y5, type="text")

model.arm.sp.y10 <- lm(log(arms_next_year_fixed+1)~terminated_y10+dcaAnyV1+US_ally+population + GDP+factor(year), data=att.2_merged)
stargazer(model.arm.sp.y10, type="text")                 

stargazer(model.arm.sp.y3, model.arm.sp.y5, model.arm.sp.y10, type="text")  # termination within 3, 5, 10 or never do not matter


# Model: DV=trade / IV=terminated_y3,5,10

model.trade.sp.y3 <- lm(log(trade_next_year+1)~terminated_y3+US_ally+population + GDP, data=att.2_merged)
stargazer(model.trade.sp.y3, type="text")  

model.trade.sp.y5 <- lm(log(trade_next_year+1)~terminated_y5+US_ally+population + GDP, data=att.2_merged)
stargazer(model.trade.sp.y5, type="text")

model.trade.sp.y10 <- lm(log(trade_next_year+1)~terminated_y10+US_ally+population + GDP, data=att.2_merged)
stargazer(model.trade.sp.y10, type="text")                 

stargazer(model.trade.sp.y3, model.trade.sp.y5, model.trade.sp.y10, type="text")

## =============================================================                  
## Reverse Causality                   
## =============================================================                  

# Check for NAs, Infs, or NaNs in your predictors
summary(att.2_merged[, c("usa_sp_times", "prev_year_arms", "US_ally", "population", "GDP")])
sapply(att.2_merged[, c("usa_sp_times", "prev_year_arms", "US_ally", "population", "GDP")], function(x) sum(!is.finite(x)))

att.2_merged_clean <- att.2_merged[is.finite(att.2_merged$prev_year_arms) & 
                                     is.finite(att.2_merged$US_ally) & 
                                     is.finite(att.2_merged$population) & 
                                     is.finite(att.2_merged$GDP) & 
                                     !is.na(att.2_merged$usa_sp_times), ]

sapply(att.2_merged_clean[, c("usa_sp_times", "prev_year_arms", "US_ally", "population", "GDP")], function(x) sum(!is.finite(x)))



# Arms tranfer => SP
model.arm.sp <- lm(usa_sp_times~log(prev_year_arms+1)+US_ally+population +GDP+factor(year)+factor(country), data=att.2_merged_clean)
stargazer(model.arm.sp, type = "text")

model.arm.sp.dca <- lm(usa_sp_times~log(prev_year_arms+1)+US_ally+population + dcaAnyV1+GDP+factor(year)+factor(country), data=att.2_merged_clean)
stargazer(model.arm.sp, model.arm.sp.dca, type = "text")

# Trade => SP
model.trade.sp <- lm(usa_sp_times~log(prev_year_tradevol)+US_ally+population + GDP+factor(year)+factor(country), data=att.2_merged_clean)
stargazer(model.trade.sp, type = "text")                   
summary(model.trade.sp)                   










